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Abstract 

We analyse the nonlinear Kuramoto-Sivashinsky equation to de- 
velop accurate discretisations modeling its dynamics on coarse grids. 
The analysis is based upon centre manifold theory so we are assured 
that the discretisation accurately models the dynamics and may be 
constructed systematically. The theory is applied after dividing the 
physical domain into small elements by introducing isolating internal 
boundaries which are later removed. Comprehensive numerical solu- 
tions and simulations show that the holistic discretisations excellently 
reproduce the steady states and the dynamics of the Kuramoto-Siva- 
shinsky equation. The Kuramoto-Sivashinsky equation is used as an 
example to show how holistic discretisation may be successfully ap- 
plied to fourth order, nonlinear, spatio-temporal dynamical systems. 
This novel centre manifold approach is holistic in the sense that it 
treats the dynamical equations as a whole, not just as the sum of 
separate terms. 



Keywords: Kuramoto-Sivashinsky equation, low-dimensional modelling, 
computational discretisations 

*Dcpartment of Mathematics and Computing. University of Southern Queensland, 
Toowoomba, Queensland 4352, Australia. 

^Department of Mathematics and Computing, University of 
Southern Queensland, Toowoomba, Queensland 4352, Australia. 
http : //www .sci .usq. edu. au/ staf f /aroberts 



1 



Contents 

AMS Subj Class: 37M99, 37L65, 65M20 

Contents 



1 Introduction 



2 Use a homotopv in the inter-element coupling! 5 



2.1 Introduce internal boundaries between elements! 6 



2.2 Centre manifold theory supports the discretisation! 8 



2.3 Approximate the shape of the centre manifold! 10 



3 Various holistic models! 12 



3.1 Some holistic discretisations! 12 



3.2 Illustration of subgrid field enhances our view! 15 



3.3 The holistic discretisations are consistent! 17 



4 Holistic models accurately give steady states! 20 



4.1 Reference accurate steady states! 21 



4.2 Holistic models are accurate on coarse grids! 22 



4.2.1 Bifurcation diagrams show success! 22 



4.2.2 Holistic models outperform centered differences! .... 25 



4.2.3 Grid refinement improves accuracy! 26 



4.3 Comparison to Galerkin approximations! 27 



4.4 Coarse grids allow large time steps! 29 



5 Holistic models are accurate for time dependent phenomenal 31 



5.1 Dynamics near the steady states are reproduced! 32 



5.2 Extend the Hopf bifurcations! 34 



5.3 Dynamics of periodic patterns without odd symmetry! 36 



6 Conclusion 40 



A Computer algebra derives the discretisation! 41 



Tony Roberts, February 1, 2008 



1 Introduction 



3 




44 



1 Introduction 



The Kuramoto-Sivashinsky equation, here 

du d A u ( du d 2 u 



dt ^ dx A ^ \ dx dx 2 



) 



0. 



(1) 



was introduced by Sivashinsky jST] as a model of instabilities on interfaces 
and flame fronts, and Kuramoto ^Hj as a model of phase turbulence in chem- 
ical oscillations. It receives considerable attention as a model of complex 
spatio-temporal dynamics [IH1 Ell El e.g.]. In the form (JTJ), with 2n pe- 
riodic boundary conditions, a is a bifurcation parameter that depends upon 
the size of the typical pattern [20]. The Kuramoto-Sivashinsky equation 
includes the mechanisms of linear negative diffusion au xx , high-order dissi- 
pation 4u xxxx , and nonlinear advection/steepening auu x . The system (JH) has 
strong dissipative dynamics arising from the fourth order dissipation. Many 
modes of this system decay rapidly because of this strong dissipation. Thus 
the dynamics are dominated by a relatively few large scale modes. We create 
and explore the macroscopic modelling of the Kuramoto-Sivashinsky dynam- 
ics using holistic discretisation as initiated by MacKenzie & Roberts |18j . 

We study the Kuramoto-Sivashinsky equation here for several reasons. 
Firstly, the pde is fourth order and therefore, following the example of Burg- 
ers' equation [53], provides a further test case for the application of the 
holistic approach to higher order dissipative PDEs. Secondly, the Kuramo- 
to-Sivashinsky equation is analogous to the Navier-Stokes equations of fluid 
dynamics. Holmes, Lumley & Berkooz ^2] argued that these analogies exist 
on two levels: in the energy source and dissipation terms of both dynamical 
systems; and in the reflection and translational symmetries of the Kuramo- 
to-Sivashinsky equation and the spanwise symmetries of the Navier-Stokes 
equations in the boundary layer. This analogy between symmetries suggests 
the Fourier series and corresponding modal interactions are comparable for 
these two problems. Thirdly, Cross & Hohenberg [Sj describes how the Kura- 
moto-Sivashinsky equation exhibits the complexities of weak turbulence or 
spatio-temporal chaos. The complex dynamics of the Kuramoto-Sivashinsky 
equation is searching test of the performance of the holistic approach to 
coarse grained modelling of dynamical systems. 

Approximate inertial manifolds and variants El El HI El e.g.] cap- 
ture the long-term low dimensional behaviour of the Kuramoto-Sivashinsky 
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equation. Most constructions of approximate inertial manifolds are based 
upon nonlinear Galerkin methods [221 1201 e.g.]. Approximate inertial 
manifolds are generally constructed by finding global eigenf unctions of the 
linear dynamics. Our approach is similar to these methods in that we project 
onto natural solutions of the pde, and performs nearly as well, see ^4.31 But 
in contrast, the holistic approach undertaken here bases analysis upon the 
local dynamics within and between finite elements and thus we contend it 
will be more useful in applications; for example, it is readily adapted to the 
modelling of a wide variety of physical boundary conditions |2Zj ■ 

Our approach is to divide the spatial domain into disjoint finite elements 
(* )2.1jl . Initially these finite elements are decoupled and so dissipation causes 
the solution to exponentially quickly become constant in each element. We 
then couple the elements together so that information is exchanged between 
elements — parameterised by a coupling parameter 7 so that 7 = 1 recovers 
the original Kuramoto-Sivashinsky equation. The coupling drives the evolu- 
tion of the field in each element. We solve the Kuramoto-Sivashinsky pde 
within each element with the coupling and hence resolve subgrid scale dynam- 
ics. Centre manifold theory (HI EH e.g.] then provides the rigorous support 
for holistic models as introduced by Roberts [22] for Burgers' equation and 
discussed in 



A low order analysis, reported in of the Kuramoto-Sivashinsky equa- 
tion ([TJ favours the discretisation 

duj 4:Uj +2 — 16uj + i + 24uj — 16iij_i + 4uj_ 2 
—Uj +2 + lQuj + i — 30ua + 16iij_i — tt,_2 



a 



12h 2 



1 n ( nl Wj+l ~ Uj-l U j+1 ~ U j-1 _ U j+2 U j+ i - Uj^Uj-! } _ 

+ a \ Uj Ah + Ah 12h )~U^ Z ) 

where the UjS are grid values spaced h apart. The first two lines of the 
holistic discretisation (J2J) shows the holistic method generates conventional 
centered finite difference approximations for the linear terms Au xxxx and au xx . 
The third line details a specific nonstandard approximation for the nonlinear 
term auu x : it is a mix of three valid approximations to uu x ; the specific 
mix is determined by the subgrid scale modelling of physical processes in the 
holistic approach, see £ 13.21 The holistic discretisation is not constructed by 
discretising the Kuramoto-Sivashinsky equation (JJ) term by term, rather the 
subgrid scale dynamics of $1} together with inter-element coupling generate 
the specific holistic discretisation (J2J). 
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The discretisation (J2J) is a low-order approximation. Centre manifold the- 
ory provides systematic refinements. Analysis to higher orders in nonlinearity 
or inter-element interaction, discussed in £0 gives further refinement to the 
discretisation. The higher order terms come from resolving more subgrid 
scale interactions. These higher order analyses lead to higher order consis- 
tency, as element size h — > , between the equivalent pdes of the holistic 
discretisations, such as (J2J), and the Kuramoto-Sivashinsky pde (see fc>.3|) . 
Such consistency is further justification for our approach. 

The bulk of this paper is then a comprehensive comparative study of the 
various models of the Kuramoto-Sivashinsky dynamics; further details are 
reported by MacKenzie [19] . A detailed numerical study of the holistic pre- 
dictions for the steady states of the Kuramoto-Sivashinsky equation is the 
focal point of Section |3J followed by an exploration of the holistic predictions 
for the time dependent phenomena of the Kuramoto-Sivashinsky equation 
in Section |S| We look at: the predicted steady states, their stability and 
compare bifurcation diagrams; the dynamics near the steady states; Hopf 
bifurcations leading to period doubling sequences; and the spatio-temporal 
patterns at relatively large nonlinearity parameter a. We find that the holis- 
tic models have excellent performance on coarse grids thus leading to sim- 
ulations that may use large time steps. The excellent performance detailed 
herein is further evidence that the holistic approach is a robust and useful 
method for discretising pdes. 



2 Use a homotopy in the inter-element 
coupling 

The construction of a discretisation is based upon breaking the spatial do- 
main into disjoint finite elements and then joining them together again. We 
control this process by a coupling parameter 7 that smoothly parametrises 
the transition between decoupled elements and fully coupled elements for 
which we recover a model for the original pde. Furthermore, we construct 
the model using solutions of the pde within each element and hence resolve 
subgrid scale dynamics. Centre manifold theory jHlESl e.g.] provides the rig- 
orous support for holistic models as introduced by Roberts for Burgers' 
equation. 
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Figure 1: An example of the ID grid with regular elements of width h. The 
jth element is centered about the grid point Xj. The vertical blue lines form 
the element boundaries, which for the jth element are located at £j±i/2 = 
(j ± l/2)h. 

2.1 Introduce internal boundaries between elements 

Establish the spatial discretisation by dividing the domain into m elements 
of equal and finite width h and introducing an equispaced grid of collocation 
points, Xj = jh, at the centre of each element, see FigureHJ 1 Express the sub- 
grid field in the jth element by u = Vj (x, t) — we solve the Kuramoto-Siva- 
shinsky pde (jlj with inter-element coupling introduced via artificial internal 
boundary conditions (ibcs). We introduce a homotopy in an inter-element 
coupling parameter 7: when 7 = the elements are effectively isolated from 
each other, providing the basis for the application of centre manifold theory; 
whereas when evaluated at 7 = 1, the elements are fully coupled together and 
hence the discretised model applies to the original pde. Since the Kuramo- 
to-Sivashinsky equation is fourth order we require four iBCs for each element 
to ensure satisfactory coupling between neighbouring elements. Here we use 
the non-local IBCs 

S x Vj(x,t) = j5v j±1/2 (x,t) atx = x j±1/2 , (3) 
6lvj(x,t) = i 2 S 3 v j±1/2 (x,t) at x = Xj ±1/2 , (4) 

which are an extension of the non-local IBCs explored by Roberts j2E] for 
Burgers' equation; a local possibility for the ibcss was explored by MacKen- 
zie jinj. These non-local IBCs involve the centered difference operators S and 5 X : 
the operator S x denotes a centered difference in x only, with step h; whereas 
the operator 5 denotes a centered difference applied to the grid index j with 
step 1; so for example, the first IBCs (jHJ) is 

Vj(x j±1 ,t) - Vj(xj,t) = j[v j+1 (x j±1 , t) - Vj(xj, t)} . (5) 

Note: the field Vj(x, t) extends analytically to at least Xj±2, to allow the appli- 
cation of the non-local IBCs (@J). The physical interpretation of these IBCs is 

principle, elements may be of unequal size. However, to simplify the analysis, herein 
all elements will be of equal width h. 
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Figure 2: Schematic diagram of the fields Vj(x, t), Vj + i(x, t) and Vj_i(x, t) for 
the non-local ibcs with 7 = 1. See the fields pass through neighbouring 
grid values Uj and it,-±i, and also Wj±2 when appropriate. 

not obvious. Firstly, when 7 = 0, (jHHU) ensures the first and third differences 
in x of the field Vj centered about the element boundaries ^±1/2 are zero. 
These isolate each element from its neighbours as there is then no coupling 
between them. In each element vj(x,t) = constant is an equilibrium. It is 
dynamically attractive provided the instability controlled by a/h 2 is not too 
large compared with the dissipation of order l/h A . This simple class of piece- 
wise constant solutions provide the basis for analysing the 7 7^ case when 
the elements are coupled together. Secondly, the non-local ibcs evaluated at 
7 = 1 requires that the field Vj(x, t), when extrapolated to Xj±i and Xj±2, is to 
equal the grid point value of the subgrid field of that element, Uj±i and Uj± 2 
respectively. See the schematic representation in Figure El of these non-local 
boundary conditions evaluated at 7 = 1. This restores sufficient continuity 
to ensure the holistic model applies to the original pde. 

The inter-element coupling parameter 7 controls the flow of information 
between neighbouring elements. We construct solutions as power series ex- 
pansions in the coupling parameter 7. 2 When 0{^j 2 ) terms are neglected in 
the holistic model, the field in the jth element involves information about the 
fields in the j ± 1 elements. Similarly, when (9 (7 s ) terms are neglected in the 
approximation, the field in the j'th element involves information about the 
fields in the j ± 1 and j ± 2 elements. Consequently, the order of 7 retained 
in the holistic model controls the stencil width of the discretisation. 

2 Such homotopies are used successfully in other numerical methods. For example, 
Liao |17| proposed a homotopy in his general boundary element method from auxiliary 
linear operators whose fundamental solutions are well known. In our application the 
homotopy is only in the ibcs. 
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Roberts argued that this particular form of the non-local ibcs ensure 
that these holistic models are consistent with any given pde to high orders 
in the grid size h as h — > . 



2.2 Centre manifold theory supports the 
discretisation 



The existence, relevance and approximation theorems [HI IH e.g.] of centre 
manifold theory apply to the Kuramoto-Sivashinsky PDE with IBCs (0- 
0J). Similar to the application to Burgers' equation by Roberts [22], the result 
here is support for a low dimensional discrete model for the Kuramoto-Siva- 
shinsky dynamics at finite grid size. 

Theoretical support is based upon the piecewise constant solutions ob- 
tained when all the elements are insulated from each other. Adjoin to the 
Kuramoto-Sivashinsky pde (0) the dynamically trivial equations for the cou- 
pling parameter 7 and the nonlinearity parameter a, 



^7 da 
dt = ~dt = 



(6) 



and consider the dynamics in the extended state space 7, a). Adjoining 

such trivial equations for parameters is commonly used to unfold bifurcations 
(HI §1-5]. In this extended space there is a subspace of fixed points with 
u = constant in each element and 7 = a = . Linearizing the pde and IBCs 
about each fixed point, u = constant + u'(x, t) , gives 



dvf_ 
~dt 



■— — r , such that Sxu'(x,t)\„-„ — 5lu'(x, t)\ =0 



The nth linear eigenmode associated with each element is 



a 



7 



, u' oc e Xnt cos 



77.7T 
~h 



X 



(7) 



for the non-local IBCs (JBH1I), where n = 0, 1, . . . and the eigenvalue A n = 
— n 4 7r 4 //i 4 . There are also the trivial modes 7 = const and a = const. 
Therefore, in a spatial domain of m elements there are m + 2 zero eigenval- 
ues: one associated with each of the m elements; and two from the trivial © • 
All other eigenvalues are negative and < — 7r 4 //i 4 . Thus, the existence theo- 
rem, see jH p. 281] or P-96], guarantees that a m + 2 dimensional centre 
manifold Ai exists for the Kuramoto-Sivashinsky pde with the trivial (JHJ) 
and ibcs (jSHU)- 
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We parametrise the m + 2 dimensional centre manifold Ai by the m + 2 
parameters 7, a and the grid values Uj. 3 Denote u as the vector of the m grid 
values. Thus for some function v the centre manifold Ai is 

u(x, t) = v(x; u, 7, a) . (8) 

Equivalently, we decompose the centre manifold field into that for each ele- 
ment: 

u = v(x;u,-y,a) = u,-y, a) Xj(x) , (9) 

3 

where the characteristic function Xj( x ) is 1 when Xj_i/2 < x < Xj+1/2 , and 
otherwise; view the centre manifold as the union of all states of the collec- 
tion of subgrid fields Vj(x; u, 7, a) over the physical domain. The correspond- 
ing amplitude condition, that the field in each element has to pass through 
its grid value, is 

Uj — v(xj] u, 7, a) . (10) 

The existence theorem j3] also asserts that on the centre manifold the grid 
values Uj evolve deterministically in time according to the system of odes 

Uj = duj/dt = gj(u,j,a) , (11) 

where gj is the restriction of the Kuramoto-Sivashinsky pde (JJJ with the 
trivial © and ibcs dHH3) to the centre manifold Ai. It is this evolution 
of the grid values that gives the holistic discretisation. 

Note that the centre manifold Ai is global in u but local in 7 and a. 
When the parameters 7 = a = the Kuramoto-Sivashinsky pde has a 
m dimensional "centre" subspace £ of fixed points with the field u being 
independently constant in each element; these are fixed points for all u. 
When the parameters 7 and a are non-zero this subspace is "bent" to the 
curved centre manifold Ai. Thus the models we construct are valid for small 
enough 7 and a, although we use them at finite 7 and a, but are formally 
valid for all |w|. Numerical solutions of the centre manifold models, such as 
those in ^3.21 indicate that parameter values as large as 7 = 1 and a = 20-50 
are indeed within the range of validity of our approach, even on relatively 
coarse grids. 

We now support the claim that the evolution of the discrete grid val- 
ues fjTTJl actually models the Kuramoto-Sivashinsky system (JTJ). The rele- 
vance theorem of centre manifolds, [U p. 282] or [H21 p. 128], guarantees that all 

3 These grid values are one choice for the measure of the field u in each element. Other 
choices are possible, but the grid values appear most convenient. 
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solutions of the Kuramoto-Sivashinsky system ((TJ with © and the ibcs 
13]), which remain in some neighbourhood of the subspace £ in (u(x),7,a) 
space are exponentially quickly attracted to the centre manifold M. and 
thence to a solution of the m discrete odes ()lljl . For our application of centre 
manifold theory to the holistic model we seek regimes where this neighbour- 
hood includes 7 = 1 and a of interest. We estimate the rate of attraction 
by the leading negative eigenvalue, here Ai = — 7r 4 //i 4 . The actual rate 
of attraction may be less due to the difference between centre manifold A4 
and centre subspace S, but Ai will be the correct order of magnitude. This 
ensures the so-called asymptotic completeness [22]: after the exponentially 
quick transients of the approach to Ai by any trajectory, the evolution of 
the discretisation (fTTj) on M. accurately models the Kuramoto-Sivashinsky 
pde (P). 

2.3 Approximate the shape of the centre manifold 

Having established that we may find a low dimensional description (fH UTTj) 
of the interacting elements that is relevant to the Kuramoto-Sivashinsky 
system ((l}, we need to construct the shape of centre manifold and the cor- 
responding evolution on the manifold. 

The approximation theorem of Carr & Muncaster p. 283] assures us 
that upon substituting the ansatz (jEHHJ into the complete system and solv- 
ing to some order of error in a and 7, then M. and the evolution thereon 
will be approximated to the same order. However, we need to evaluate the 
approximations at the coupling parameter 7 = 1 because it is only then that 
the artificial internal boundaries are removed. Thus the actual error of the 
model due to the evaluation at 7 = 1 is not estimated. However, the holis- 
tic method for discretising the Kuramoto-Sivashinsky equation is supported 
three ways: firstly, the smooth homotopy from 7 = with large spectral gap 
to the gravest decaying mode with decay rate ~ — 7r 4 //i 4 ; secondly the holistic 
models are consistent with the Kuramoto-Sivashinsky PDEto high order in 
grid size h, see thirdly, we see in Sections the holistic models model 
accurately both steady state solutions and time dependent phenomena of the 
Kuramoto-Sivashinsky system. 

To construct the centre manifold, we solve for the field v j in each element. 
For definiteness, here we consider domains periodic in space, or equivalently 
elements far from the influence of any physical boundary. By translational 
symmetry of the Kuramoto-Sivashinsky pde the subgrid field in each 
element is identical, except for the appropriate shift in the grid index j. 
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Thus, we construct the subgrid field and evolution for a general jth element, 
see some examples in Sectional 

The algebraic details of the derivation of the centre manifold model (jS HTTj) 
are handled by computer algebra. In an algorithm introduced by Roberts [2*4] . 
iteration drives to zero the residuals of the governing pde (0) and its ibcs (JSJ- 
EJ) and amplitude condition (|T0|) . Since the algebraic details of the construc- 
tion are tedious, they are not given; instead see the computer algebra proce- 
dure of |2S]. 

This computer algebra is based upon driving the residuals of the governing 
equations to zero in the following manner. Recall from §2.21 that the centre 
manifold (JHJ) is parametrised by the grid values u and that the evolution of 
the grid values is given by (fTTj). Thus substitute these into the Kuramoto- 
Sivashinsky pde and seek to solve 



8va v-^ dvj , d^Vj ( d 2 VA dvi , 

3 -9k = - <* —2 + ^.-2 , (12) 



dt ^-^ duk- yk dx 4 " I dx 2 ^ 1 dx 



- J-k 

k 



together with the non-local ibcs (j^lEJ) and the amplitude equation (|10[). to 
some order in parameters 7 and a. The iteration is that given any approx- 
imation, denoted by ~ we seek corrections, denoted by primes, such that 
Vj = Vj + v'j and gj = gj + g'j, better satisfy the Kuramoto-Sivashinsky pde. 
Thus in each iteration we solve a problem of the form, 



-4—^ = o' + Residual , (13) 
ox 4 J 

„ , , v-^ 9va j d 4 v; ( d 2 Vj „ dv- \ . . 

Residual = ^-± gt+i -i +a ^-i + Vj -i j , (14 ) 

together with the ibcs, for the corrections, primed quantities, to the subgrid 
field and the evolution of the grid values. Note: the residual in |H}) is the 
residual of the Kuramoto-Sivashinsky system for the current approximation. 
The iteration scheme starts with the linear solution in each element, namely 
Vj(x,u,'j,a) = Uj and gj(u, 7, a) = 0. The iteration terminates when the 
residuals of the Kuramoto-Sivashinsky pde |T2|) . and the ibcs, are zero to 
some order in 7 and a. Then theory assures us that the subgrid field in each 
element and the evolution of the grid values are correct to the same order in 
7 and a. 



where 
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3 Various holistic models 



Here we record holistic models of the Kuramoto-Sivashinsky pde ((TJ), to 
various orders in coupling parameter 7, governing the width of the numerical 
stencil, and in the nonlinearity parameter a. For use the models need to 
be evaluated at 7 = 1 as then the non-local ibcs (JBH3J ensure sufficient 
continuity in the solution field. We write the models in terms of the centered 
difference and mean operators, 

5uj = Uj+x/2 - Uj-1/2 and fiuj = (ttj+1/2 + u j _ 1 / 2 )/2 , 

respectively. The models are constructed using a reduce program adapted 
from [2B]- We only present in detail here holistic models to errors (9 (a 2 ) as 
the level of complexity increases enormously with the order of a. 



3.1 Some holistic discretisations 



In order to represent the spatial fourth derivative in the Kuramoto-Sivashin- 
sky equation, we need at least a 5 point stencil approximation. Thus we de- 
termine the interactions between at least next-nearest neighbouring elements 
by obtaining up to at least quadratic terms in the coupling parameter 7. 



The C?(7 3 ,a: 2 ) holistic discretisation is 

7«c 2 7« x 47 2 4 7 2 "e 4 
Uj = —-rrrO U,: —UjOUUj — Uj H — Uj 

3 j h h I2h 2 3 

+ — (2uj5 3 /iUj + 5 2 Uj5 3 fiUj + 5 4 Uj5fiUj) 

+ 0( 7 3 ,a 2 ), (15) 

for the non-local ibcs (JBJ-HJ). This forms a basic 5 point stencil approximation, 
since the evolution itj involves just Uj, Uj±i and Uj± 2 . The first line of (fT3j) . 
when evaluated at 7 = 1 , gives a 2nd order centered difference approximation 
for the hyper diffusion term 4:U XXXX , a 4th order centered difference approxi- 
mation to the linear growth term au xx , and a 2nd order centered difference 
approximation to the nonlinear advection term auu x . The second line modi- 
fies the nonlinear discretisation to account for interaction with effects caused 
by the next-nearest neighbour elements. 

The holistic discretisation (fTK|) contains the approximation 

.... I _ / Mj+l ~ tf+l - U j+2 U j+1 - Mj- 2 Mj-A , , 

UUx ^ ~ \ Uj 4h + 4h 12h J ■ [ b) 
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when evaluated at 7 = 1. This is a 1/2 : 1 : — 1/2 mix of the approximations 

.... I ~ .. %+l ~ ~ ~ ~ ^+2^+1 ~ ^-2^-1 n? s 

respectively. This particular nonstandard approximation (llfij) to the non- 
linear term auu x , arises due to the modelling of subgrid scale interactions 
between the Kuramoto-Sivashinsky equation and the inter-element coupling. 
Such nonstandard approximations generated through this approach can have 
robust numerical characteristics [2E1 §2]. 

The 0(7 4 ,ai 2 ) holistic discretisation is 

7« r2 l a r 4 7 2 ,4 , 7^ r4 
It 7 = — TTrO Mi —UjOUUj — Uj H — rO Wj 

J /i 2 J h 3 p J /i 4 3 12/i 2 J 

27 3 rfi 7 3 a efi 

3/i 4 J 90/i 2 3 

7^ 'o? 

+ — (2uj5 3 fiUj + 5 2 Uj5 3 fiUj + 5 4 Uj5uuj^ 



7 a 



(16 u/ 5 + 30 8 4 UjS 3 fj,Uj + 40 S 2 u j S 3 fiu j 



480/i 

+ 40 5 UjSfiUj + 28 S 2 Uj5^fiUj + 14 5 6 Uj5fiUj 

+ 76 4 Uj8 5 fMj + iPu^nuj) + C ( 7 4 , a 2 ) . (18) 

This discretisation forms a 7 point stencil approximation, involving itj, Uy±i, 
Mj±2 and Wy±3- The first two lines of IfTSjl . when evaluated at 7 = 1, give 
a 4th order centered difference approximation to the hyper diffusion term, a 
6th order centered difference approximation to the linear growth term, and 
a 2nd order centered difference approximation to the nonlinear advection 
term. The third and remaining lines account for higher order subgrid scale 
dynamics of the nonlinearity and its inter-element coupling to generate a 
4th order centered difference approximation to the nonlinearity uu x . 

The C?(7 5 ,a: 2 ) holistic discretisation is 

7«, 2 7« r 4 7 2 . 7 2 a 4 
Uj = -—5 Uj - —UjdfiUj - —6 Uj + Uj 

2 7 3 6 7 3 a 6 

+ W 6 U3 ~ 9W 5 

7 7 4 rs„. , 7 4 « X 8 



6O/1 4 3 560/i 2 
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+ — (2uj5 3 fiUj + 5 2 Uj5 3 fiUj + 5 4 Uj5fiUj) 

7 

— (16 Uj5 5 fiUj + 30 5 4 Uj5 3 fiUj + 40 5 2 Uj5 3 fiUj 

+ 40 5 A Uj5{iUj + 28 5 2 u j 5 5 j2u j + 14 5 6 Uj5fj,Uj 
+ 75 4 Uj5 5 fiUj + 7S 6 Uj5 3 /iuA 

4 

(432u/ 7 /i% + 3528 5 2 u j 6 5 fm j + 1507 5 2 Uj5 7 ^ 



60480/i 

+ 3780 S 4 Uj5 3 fiUj + 3951 S 4 Uj5 5 fiUj + 984 S A Uj^ 

+ 1764 5 6 w i 5/iw i + 3419<5 6 M j 5 3 /i^ + 1414 5 6 Mj 5 5 /i% 

+ 164 5 & u j 5 7 ^u, J + 523 5 8 u j 5fiu j + 656 5 s u j 5 3 fiu j 

+ 164 5*Uj5*nu 3 ) + (7 5 , a 2 ) . (19) 

This forms a 9 point stencil approximation, involving only it.,-, Uj±i, %±2 ; 

and Mj±4- The first two lines of (fT9"j) when evaluated at 7 = 1 give a 
6th order centered difference approximation for the hyperdiffusion term, an 
8th order centered difference approximation for the linear growth term and 
a 2nd order centered difference approximation for the nonlinear advection 
term. The third and remaining lines provide modifications to model the 
nonlinear uu x to 6th order through resolving subgrid scale dynamics. 

We do not code these discretisations manually. Instead, the computer al- 
gebra program at [2Hj is used with the UNIX editor sed to automatically write 
the discretisation in a form suitable to be input to MATLAB for numerical 
exploration. 



Compare to conventional centered difference models. Traditional 
direct finite differences generate the following approximations to the Kura- 
moto-Sivashinsky pde ((TJ: 

• 5 point, 

= —^u j Sfiu j -—Su j -—Su j + 0(h 2 ); (20) 



7 point, 



a ( r 1 r3 A a (si 1 U 
Uj = — — UjOUUj U~0 UUj ) — —r I Uj Uj 

3 H J P 3 6 3 P 3 J h 2 V 3 12 3 



4 

h 1 



(s 4 Uj - jU 6 «jJ +c (^ 4 ) ; (21) 



Tony Roberts, February 1, 2008 



3 Various holistic models 



15 



• 9 point, 




Consider the different view of the errors for the discretisations: the cen- 
tered difference approximations ([20H22)) are justified by consistency as grid 
size h — > ; whereas the holistic discretisations (fTBTfT^j) are supported by 
centre manifold theory at finite grid size h. The errors in the centre mani- 
fold approach are due to the truncation of dependence in the inter-element 
coupling parameter 7 and the nonlinearity parameter a. However, as argued 
by Roberts [23 for linear systems and as demonstrated in H'd.'Sl the particu- 
lar choice of the ibcs (jHHU) ensures that the holistic discretisations are also 
consistent as h — > with the Kuramoto-Sivashinsky pde ((H). 

3.2 Illustration of subgrid field enhances our view 

Recall the collection of subgrid fields Q over the physical domain form a 
state on the centre manifold. Here we plot some example subgrid fields 
for various holistic models. In particular, we examine subgrid fields of the 
holistic models of steady states of the Kuramoto-Sivashinsky pde Q at 
nonlinear parameter a = 20 and a = 50 . This is intended to reinforce the 
link between the abstract centre manifold description of the dynamics and 
the physical subgrid fields for the low order holistic models. We compare 
the fields with the Lagrangian interpolation that underlies traditional finite 
differences. Recall that the key methodology difference is that the subgrid 
fields of the holistic models are constructed by actual solutions of the Kura- 
moto-Sivashinsky pde, see < j2.3l 

We restrict attention to odd symmetric solutions that are 27r-periodic. 
This is done to compare with the numerical investigations of Jolly ^1] which 
we consider in more detail in Sections 0] and 03 Set a grid of 8 equi-spaced 
elements on the interval [0, ir]. The subgrid fields are plotted for approxi- 
mations to the steady states of the Kuramoto-Sivashinsky equation (pj with 
these periodic boundary conditions, computed using holistic discretisations 
at a = 20 and a = 50 . 
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Figure 3: Subgrid field (green curve) of the holistic model flT5)) and a La- 
grangian interpolant (magenta curve) constructed through a 2nd order cen- 
tered difference approximation for a steady state of the Kuramoto-Sivashin- 
sky equation at a = 20, with 8 elements on [0,7r]. An accurate solution is 
also plotted in blue. 

Figure 4: Subgrid fields of the holistic models with errors 
(green), 0(7 4 , a 2 ) (fTHj) (olive green) and 0(j 5 ,a 2 ) (fTTij) (cyan), for a steady 
state of the Kuramoto-Sivashinsky equation at a = 20, with 8 elements on 
[0,7r]. An accurate solution is also plotted in blue. 

Figure El displays an accurate solution (blue curve) of the Kuramoto- 
Sivashinsky pde to compare with the subgrid field (green curve) of the 
5 point stencil O (7 s , a 2 ) holistic approximation ()15j) (green discs), and the 
Lagrangian interpolant (magenta curve) constructed through a 2nd order cen- 
tered difference approximation (magenta discs), for a steady state at a = 20 . 
Observe the collection of subgrid fields forms the field u which is a state on 
the centre manifold. The subgrid field of the holistic model more accurately 
represents the steady state of the Kuramoto-Sivashinsky equation at a = 20 , 
on this coarse grid. 

Higher order holistic models improve the accuracy and continuity of the 
subgrid field. Figure |U displays the subgrid fields of three holistic models for 
the same steady state of the Kuramoto-Sivashinsky pde depicted in Figure El 
for a = 20. The (9 (7 s , a 2 ) holistic model ()15|) (green) is the least accurate 
and has the largest jump at element boundaries. The (9(7 4 ,a 2 ) (fT8l) model 
(olive green) displays improvement over the holistic O (7 s , a 2 ) approxima- 
tion. The (D(y 5 ,a 2 ) (TTTjj) model (cyan) is the most accurate, being almost 
indistinguishable from the correct curve. 

Figure El shows a steady state of the Kuramoto-Sivashinsky pde at a = 
50 . The accurate field is symmetric (blue curve). For this value of the nonlin- 
earity there is no steady state solution for centered difference approximations 
of either 2nd (|2T)|) . 4th (f2~Tj) or 6th order (|2*2*|) on this coarse grid of 8 elements 
on [0, 7r]. However, the 5 point stencil holistic approximation with errors 

Figure 5: Subgrid fields of the holistic models with errors 0(7 3 ,a 2 ) ATS]) 
(green), C(7 4 ,a 2 ) (JTHJ) (olive green) and 0(^j 5 ,a 2 ) (JEJ) (cyan), for a steady 
state of the Kuramoto-Sivashinsky equation at a = 50, with 8 elements on 
[0,7r]. An accurate solution is also plotted in blue. 
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(9 (7 s , a 2 ) ()15j) (green) models this steady state of the Kuramoto-Sivashin- 
sky equation even for such a large value of the nonlinearity on this coarse 
grid. This O (7 s , a 2 ) holistic solution has significant jumps across the sub- 
grid field at element boundaries; moreover, the subgrid field is not symmetric 
and is most inaccurate near the centre of the spatial domain considered here. 
The 7 point stencil holistic approximation with errors (9(7 4 ,a 2 ) (|18|) (olive 
green) is more accurate with smaller jumps between neighbouring the subgrid 
fields, but is also not symmetric. The 9 point stencil holistic approximation 
with errors 0(^j 5 , a 2 ) (fT9*|) (cyan) is the most accurate of the holistic models 
illustrated here; it is symmetric and the jumps between neighbouring subgrid 
fields are almost indiscernible. 

These illustrations of the subgrid fields of steady states of the Kuramo- 
to-Sivashinsky equation at a = 20 and a = 50 indicate the holistic models 
perform well even at such large values of a supposedly small parameter. The 
performance of the holistic models are explored further in Section|3]for steady 
states and Section El for time dependent phenomena. 

3.3 The holistic discretisations are consistent 

Holistic models constructed by implementing the ibcs (j^HU have dual justi- 
fication : they are supported by centre manifold theory for small enough 
a and 7; as well as being justified by their consistency as the grid size h 0. 
We explore consistency as a well established feature of numerical analysis. 4 

Here we examine the equivalent pdes for the holistic discretisations (fToT - 
fl"9j) evaluated at 7 = 1, and the centered difference approximations (|20Tl22|) . 
These equivalent pdes establish the (D(h 2p ~ 2 ) consistency with the Kuramo- 
to-Sivashinsky pde for holistic models constructed with residuals (D(y p+1 ). 

Roberts j2Hj proved that using ibcs of the form introduced in £|2] and 
retaining terms up to 7 P in the holistic approximations results in approxima- 
tions which are consistent with the linear terms of the Kuramoto-Sivashinsky 
equation (JTJ) to C(/i 2p ~ 2 ), provided p > 2. However, it appears that using 
the ibcs dHHU) also ensures 0(h 2p ~ 2 } consistency for the nonlinear terms. As 
yet no formal proof exists of this nonlinear consistency, but all holistic mod- 
els of the Kuramoto-Sivashinsky equation, containing terms up to 7 7 and a 4 
and constructed using (J3HIJ are nonlinearly consistent (although not all are 
recorded here). 

4 But note that high order consistency is not a primary goal of this holistic approach, 
since we aim to develop and support models for finite element size h. 
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Find the equivalent pdes for the various discretisations by expanding the 
discretisations in grid size h about a grid point Xj. That is, write 

, duj 9 h 2 d 2 Uj ^ h h k d k U4 , , 

u i±m = u, ± mh^ + » ! T ^ + E(±")* ¥ , (23) 

fc=3 

to whatever order in h is required. Computer algebra performs the tedious 
details. 



The equivalent PDE for the 5 point holistic discretisation (fTHjl . 

which retain terms up to 7 2 , is 

du ( du d 2 u\ d 4 u 2h 2 d 6 u h 4 d 8 u 

-a u— -\- — — ) 4— - — — - + 



dt \ dx dx 2 J dx 4 3 dx 6 20 dx 8 

1 d 3 ud 2 u 1 d 4 udu 1 <9 5 w 1 d 6 u 



,48<9x 3 &r 2 48<9x 4 <9x 30 ftr 5 90 dx 6 , 
+ 0(/i 6 ) . (24) 

The equivalent PDE for the 5 point centered difference approximation (|20|) is 

du ( du d 2 u\ d 2 u 

-a [ u— + — ) -4- 



h 2 I a-— -— — h a——— + -- 



l<9 3 w<9w 1 d 4 u 2d 6 



6 cte 3 dx 12 <9x 4 3 <9ic 6 
+ 0(h 4 ) . (25) 

Observe that both equivalent pdes f!24ll23|) are 0(h 2 } accurate. The coeffi- 
cients of the error terms are different in both of the these equivalent pdes, 
with those of (fl3|) having fewer error terms. 

The equivalent PDE for the 7 point holistic discretisation IjIBj) . 

which retains terms up to 7 3 , is 

du ( du d 2 u\ d 4 u 7h 4 d s u 13h 6 d w u 

— = —a I u 1 — 4 1 

dt \ dx dx 2 J dx 4 60 dx 8 756 dx 10 

,fYI_d^u^u l_^u_&\i_ 3 d 6 u du 

~ a V640 dx 4 dx 3 + 384 dx 5 dx 2 + 320 dx 6 dx 

1 d 7 u 1 d 8 u 

-I u 1 

140 dx 7 560 9x 8 
+ 0(h 8 ) . (26) 
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Whereas the equivalent pde for the 7 point centered difference approxima- 
tion (jni is 

du ( du d 2 u\ ^) 2 u 

dt \ dx dx 2 J dx 2 

4 / 1 d 5 udu 1 d 6 u 7 d 8 u\ 

\ a wl^lh +a wlw + mdx^) 

+ 0(h 6 ) . (27) 

The two equivalent pdes (I2fihl2~7j) are (9(/i 4 ) accurate, and again the holistic 
discretisation has fewer errors. 



The equivalent PDEs for the 9 point holistic discretisation (fTHj) . 

which retain terms up to 7 4 , is 

du ( du d 2 u\ d*u Alh 6 d 10 u 13h 8 d 12 u 

-a u— + — - - 4- 



dt V dx dx 2 J dx A 1890 dx 10 2700 dx 12 

8 / 3433 d 5 u d 4 u 5927 d 6 u d 3 u 499 d 7 u d 2 u 
+ a V 138240 dx 5 <9x 4 + 322560 dx 6 dx 3 + 53760 dx 7 dx 2 

29 d 8 udu 1 d 9 u 1 d 10 u\ 
+ 8^dx 8 dx + 630 U dx^ + 3150^° J 

+ 0{h 10 ). (28) 
The equivalent pde for the 9pt centered difference approximation (|2~2"J) is 

du ( du d 2 u\ d 2 u 

-a u— + -— - 4- 



dt \ dx dx 2 J dx 2 

d 7 udu 1 d 8 u 41 d w u 
— h a—-r-^;-^ — h a 



140 dx 7 dx 560 9x 8 1890 dx 10 , 
+ C(/i 8 ) . (29) 

Again the the equivalent PDEs (j28fJ2T?|) are (9(/i 6 ) accurate, with the holistic 
discretisation having fewer error terms. 

Although there is no proof of nonlinear consistency in general, we have 
demonstrated it here for these three holistic discretisations, and have found 
nonlinear consistency for all models investigated. 
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Figure 6: Accurate bifurcation diagram < a < 70 for the Kuramoto-Siva- 
shinsky equation, using a 6th order centered difference approximation with 
48 points on the interval [0, ir]. A signed L 2 norm is plotted against a 

4 Holistic models accurately give steady 
states 

The relevance of our holistic models is rigorously supported by centre mani- 
fold theory for sufficiently small parameters 7 and a . However, the holistic 
models must be evaluated at coupling parameter 7 = 1 to model the dynam- 
ics of the Kuramoto-Sivashinsky equation. The important question: Does 
evaluating the holistic models at 7 = 1 provide useful and accurate numerical 
models? Numerical experiments detailed in this and the next section provide 
strong support that it does. 

In this section we explore the accuracy of the holistic models by construct- 
ing and comparing bifurcation diagrams of the various holistic discretisations 
to conventional explicit centered difference approximations and to the bifur- 
cation diagrams presented by Jolly et al. ^3] for various traditional Galerkin 
and nonlinear Galerkin approximations. 

We restrict exploration to solutions that are both 2ir periodic and odd: 
thus 

u(x, t) = u(x + 2n, t) and u(x, t) = — u(2n — x, t) . (30) 

We also restrict the nonlinearity parameter to the range < a < 70 . These 
restrictions are to compare our results to those of Jolly et al. ^3] for approx- 
imate inertial manifold methods. For this range of nonlinearity a the trivial 
solution u = undergoes pitchfork bifurcations at a = 4, 16, 36, 64 leading 
to the unimodal, bimodal, trimodal and quadrimodal branches respectively, 
see the bifurcation diagram Figure |H1 

Such bifurcation diagrams usefully summarise qualitative and quantita- 
tive information for a large range of the nonlinearity parameter a. We use 
the software package xppaut jH], which incorporates the continuation soft- 
ware auto [7j, to calculate the bifurcation information. The information is 
then filtered through a function written in MATLAB to draw the bifurcation 
diagram. The input to xppaut is a text . ode file describing the set of odes. 
Because the holistic models contain a large number of terms the . ode files are 
generated automatically using reduce and matlab which also incorporates 
the odd periodic requirement (}3T)j) . see [THj for more details. 
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Figure 7: Some examples of the stable equilibria of the Kuramoto-Sivashin- 
sky equation. Dark blue curves are solutions along the negative unimodal and 
bimodal branches. Light blue curves are stable solutions along the negative 
trimodal branch. 

4.1 Reference accurate steady states 

Here we introduce accurate solutions for the steady states of the Kuramo- 
to-Sivashinsky equation ([TJ over the range < a < 70 as summarised in 
the bifurcation diagram of Figure HI Accurate solutions are produced by 
a 6th order accurate centered difference approximation (|22|) with 48 grid 
points on the spatial interval [0, 7r]. These provide the reference for the 
approximations on coarse grids, and serve to also introduce the conventions 
we adopt in bifurcation diagrams. 

For all the bifurcation diagrams a signed solution norm is plotted against 
the nonlinearity parameter a. This is different to the convention adopted by 
Jolly et al. [Hj but empowers us to investigate more detail by showing positive 
and negative branches — stability differs along these branches. For example, 
see in Figure El that the negative bimodal branch is stable for 16.140 < a < 
22.556, whereas the positive bimodal branch is unstable. The solution norm 
is signed corresponding to the sign of the the first grid value, u\ = u(xi) . The 
blue curves are branches of stable fixed points and the red curves are branches 
of unstable fixed points. The open squares denote pitchfork bifurcations and 
the black squares denote Hopf bifurcations. 

The labeling scheme used in Figure El follows that of Jolly et al. [TJ] and 
Scovel [3U] with the addition of a plus or minus sign depending upon the 
sign of u±. For example, the secondary bifurcation on the negative bimodal 
branch is labeled i?2^i — from the labeling scheme of Scovel with the addition 
of the — sign because it occurs on the negative branch. Figure El appears 
to show several discontinuities. For example, the positive unimodal branch 
ends at approximately a = 12 . This apparent discontinuity arises due to the 
convention adopted here of taking the sign of U\ to sign the norm: actually 
there is a continuous transformation as the positive unimodal branch and 
the negative unimodal branch transform into the negative bimodal branch. 
It is straightforward to sign the branch near the trivial solution, but away 
from the trivial solution the distinction between positive and negative may 
be ambiguous and occasionally leads to jumps in the bifurcation diagram. 

For later comparison see in Figure [7| some of the stable equilibria of the 
Kuramoto-Sivashinsky equation in the regime of interest, < a < 70 . Fig- 
ures HK,b,c show solutions on the negative unimodal branch at a = 1,5, 10 
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Figure 8: Some accurate solutions plotted with holistic and centered differ- 
ence approximations on coarse grids. Blue curves are accurate solutions, 
green curves are the holistic approximation with ibcs (JBHIJ) with errors 
(9(7 5 ,a 2 ) on 8 elements. Magenta curves are a 6th order centered differ- 
ence approximation with 8 grid points. 

Figure 9: Bifurcation diagrams for coarse grid approximations with 8 ele- 
ments on [0, 7r] for (a) holistic model 0(7 5 , a 2 ), (b) centered difference 6th or- 
der. 

respectively. Figures[7|i,e,f show solutions on the negative bimodal branch at 
a = 20, 30, 40 respectively. The dark blue curves in Figures [7g,h,i show solu- 
tions on the negative bimodal branch and the light blue curves are solutions 
on the negative trimodal branch at a = 50, 55, 60 respectively. 

4.2 Holistic models are accurate on coarse grids 

We begin investigating the performance of the holistic models by considering 
the (9(7 5 ,a 2 ) holistic model (|T9"j) (9 point stencil, (9(/i 6 ) consistent). We 
investigate its reproduction of the steady states of the Kuramoto-Sivashin- 
sky system using coarse grids on the interval [0, ir}. 

Figure |H1 shows accurate solutions of the Kuramoto-Sivashinsky equa- 
tion (P) with odd boundary conditions (jHOj) in blue. The holistic model with 
errors 0(7 5 , a 2 ) and 8 elements is shown in green: Figure |H|;,h,i, the bottom 
row, shows that the holistic model with errors 0(7 5 , a 2 ) gives at large non- 
linearity the stable bimodal and trimodal solutions, a = 50 and a = 55 , and 
the stable bimodal solution at a = 60 . Magenta curves are solutions of the 
6th order centered difference approximation (|2~2~|) with 8 grid points — it has 
equal stencil width to the holistic model. The 6th order centered difference 
approximation does not give any stable solutions for a > 50 . The holistic 
model provides reasonable solutions where comparable traditional methods 
do not. 

4.2.1 Bifurcation diagrams show success 

Now turn to the bifurcation diagram to obtain a more comprehensive view. 
We see the holistic model has good bifurcation diagrams on a coarse grid of 
8 elements and even with just 6 elements. 
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Figure 10: Bifurcation diagrams for coarse grid approximations with 6 ele- 
ments on [0, 7r] for (a) holistic model 0(7 5 , a 2 ), (b) centered difference 6th or- 
der. 

Figure 03 shows a side by side comparison of the holistic model with er- 
rors 0(7 5 , a 2 ) with 8 elements on [0, ir] and the 6th order centered difference 
approximation with 8 grid points on [0,7r]. These approximations are both 
9 point stencil approximations. The accurate bifurcation diagram is also 
plotted in grey but without any stability information. The signed L 2 norms 
for the bifurcation diagrams on the coarse grid of 8 elements are adjusted 
by a factor of \/6 to allow comparison to the accurate bifurcation diagram 
constructed with 48 grid points on [0, ir]. Throughout this paper when com- 
paring bifurcation diagrams of different grid resolutions, the signed L 2 norms 
are adjusted this way to provide a consistent reference. Figure shows the 
(9 (7 s , a 2 ) holistic model gives good agreement with the accurate bifurcation 
diagram for a < 40 and qualitatively reproduces most of the bifurcation 
picture for 40 < a < 70 . The O (7 s , a 2 ) holistic model does not detect 
the bifurcation points R^^. on this coarse grid and the bifurcation points 
i?3ti± are incorrectly identified as fold points. However, the 0(7 5 , a 2 ) holis- 
tic model finds all of the other bifurcation points in this range of a. Figure^ 
shows the 6th order centered difference approximation gives good agreement 
with the accurate bifurcation diagram only for a < 20 and qualitatively re- 
produces the bifurcation diagram for 20 < a < 40 . The 6th order centered 
difference approximation performs poorly for a > 40 . Tabled lists the values 
of a at which the bifurcation points occur and confirms the 0(7 5 , a 2 ) holis- 
tic model performs more accurately than the 6th order centered difference 
approximation on this coarse grid of 8 elements. 

Figure El is a side by side comparison of the same £>(7 5 ,a 2 ) holistic 
model to the 6th order centered difference approximation, on an even coarser 
grid of just 6 elements. The superior performance of the holistic model is 
again evident. We conjecture that the superior performance of the holistic 
discretisation is due to its systematic modelling of the subgrid scale processes. 
These bifurcation diagrams, Figures 191 and ITUl give excellent support to the 
holistic approach to generating approximations for the Kuramoto-Sivashin- 
sky equation. 

We also investigate various holistic models for the Kuramoto-Sivashinsky 
PDE by comparing bifurcation diagrams of holistic models of higher orders. 
We examine bifurcation diagrams for holistic models with errors (9(7 p ,a 9 ), 
for 3 < p < 5 and 2 < q < 4 , and find that retaining terms of higher order in 
coupling parameter 7, corresponding to wider stencil approximations, gives 
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Table 1: a values at which bifurcation points occur for the various coarse 
grid approximations; * denotes bifurcation point identified as fold point. 
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Centered 12 pts 



2nd order 15.77 21.68 48.33 57.63 34.36 44.70 58.34 58.36 

4th order 16.12 22.37 51.74 62.33 35.98 48.62 63.11 62.98 



much greater improvement in accuracy than retaining terms of higher order 
in the nonlinearity parameter a. 

Figure E] shows the bifurcation diagrams for the holistic models up to 
and including the £?(7 5 , a 4 ) holistic model. Surveying across the columns 
of Figure ^2 see the bifurcation diagrams for holistic models of increasing 
order of coupling parameter 7, corresponding to approximations of increasing 
stencil width. For example, Figure ITTk.b.c shows the bifurcation diagrams 
for the holistic models ()18j1 and (J 19)1 respectively. Surveying down the 

rows of Figure ^2 see the bifurcation diagrams for increasing orders of the 
nonlinearity parameter a. Figure ^2 illustrates the improvement in accuracy 

Figure 11: Bifurcation diagrams for the holistic models with 8 elements on 
the interval [0, n] up to and including the (9(7 5 , a 4 ) holistic model. 
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of the higher order holistic models. Note first the dramatic improvement in 
accuracy gained by moving from left to right across Figure ^JJ corresponding 
to approximations of higher orders in the coupling parameter 7. 

Second, see that less improvement is gained by moving from top to bot- 
tom of Figure El corresponding to approximations of higher order in the 
nonlinearity parameter a. There are some peculiarities about this series of 
bifurcation pictures for holistic models of increasing order in a. For the 
5 point stencil approximations displayed in the first column of Figures El 
higher orders in a appear to gain some improvement. In particular Fig- 
ures ITTH.g show the C(7 3 , a 3 ) and C(7 3 , a 4 ) holistic models reproduce the 
unstable trimodal branches that were missing from the 0(7 3 , a 2 ) bifurcation 
diagram shown in Figure fTTk. However, for the 7 point stencil approxima- 
tions displayed in the second column of Figure ^TJ holistic models of higher 
orders in a lose some features of the Kuramoto-Sivashinsky system. The 
correct behaviour of the unstable trimodal and quadrimodal branches is re- 
produced for the (9(7 4 ,a 2 ) model shown in Figure ITTb. but not reproduced 
for the higher order (9(7 4 ,a 3 ) and (9(7 4 ,a 4 ) models shown in Figures ITTfe.h 
respectively. For the 9 point stencil approximations, displayed in the third 
column of Figures ITTl the (9(7 5 ,a 2 ) holistic model shown in Figure ITTb. re- 
produces the unstable trimodal branch whereas the higher order 0[y 5 ,a 3 ) 
model shown in Figure ITTF. does not reproduce the unstable trimodal branch. 
These peculiarities suggest that while we have observed excellent performance 
of the holistic models constructed with the non-local ibcs on coarse grids, it 
may be possible that modifications could be made to the non-local ibcs such 
that higher order approximations in the nonlinear parameter are improved. 
Exploration of possible such modifications are left for further research. 

4.2.2 Holistic models outperform centered differences 

In §4.2.11 we saw that the performance of the 0{^j 5 ,a 2 ) holistic model |T9|) 
constructed with non-local ibcs was far superior to the explicit 6th order cen- 
tered difference approximation (j22j) . To complete the comparison of holistic 
models to explicit centered difference schemes, we compare the (7 s , a 2 ) JSJ) 
and 0(7 4 , a 2 ) (HHJ) holistic models to the 2nd order (jUJ) and 4th order (|2"Tjl 
centered difference approximations respectively; these are 5 point and 7 point 
discretisations respectively. 

The first row of Figure IT21 is a side by side comparison of the O (7 s , a 2 ) 
holistic model and the 2nd order centered difference approximation with 8 el- 
ements on [0, 7r]. The second row of Figure IT2*1 is a side by side comparison of 



Tony Roberts, February 1, 2008 



4 Holistic models accurately give steady states 



26 



Figure 12: Bifurcation diagrams for (a) (9(7 3 ,a 2 ) holistic model, (b) 2nd or- 
der centered difference, (c) 0(7 4 , a 2 ) holistic model and (d) 4th order cen- 
tered difference all with 8 elements on the interval [0, 7r] 

Figure 13: Bifurcation diagrams for the holistic models with 12 elements on 
the interval [0,7r]. Compare with Figure HP with 8 elements. 

the C?(7 4 , a 2 ) holistic model and the 4th order centered difference approxi- 
mation on the same coarse grid. The accurate bifurcation diagram is plotted 
in grey without any stability information. 

Although comparing Figures IT2*b.d shows some improvement is gained by 
taking higher order centered difference approximations, this improvement is 
not as pronounced as for the holistic models on this coarse grid as shown in 
Figures IT2*k.c. Both the 2nd order and 4th order centered difference approxi- 
mations fail to reproduce the correct behaviour of the unstable trimodal and 
quadrimodal branches. In contrast, even the 5 point stencil (9(7 3 ,a 2 ) holis- 
tic approximation qualitatively reproduces the trimodal and quadrimodal 
branches on the same coarse grid. The values at which the bifurcation points 
occur are listed in Table Q and confirm these holistic models outperform the 
centered difference approximations on this coarse grid of 8 elements on [0, tc]. 

4.2.3 Grid refinement improves accuracy 

Since the equivalent pde's, (|2*4^). (|2*Bj) and (|2*Hj) . for our holistic models are 
of 0[h 2 ), C(/i 4 ) and 0(/i 6 ) respectively, grid refinement should result in 
improved accuracy. 

Figure EH shows the bifurcation diagrams of the holistic models up to 
and including the 0{y 4 ,a 3 ) model on a finer grid of 12 elements on [0,7r]. 
Compare Figure [TBI with Figure ^Tj to confirm the improved accuracy for the 
holistic models on this refined grid. Tabled also shows the bifurcation points 
are more accurately reproduced for the holistic models on this refined grid. 

Figure E] is a side by side comparison of the bifurcation diagrams of 
the (9(7 4 ,a 2 ) holistic model and the 4th order centered difference approx- 
imation (j2TJ)- The accurate bifurcation diagram is shown in grey. See the 

Figure 14: Bifurcation diagrams for (a) (9(7 4 ,a 2 ) holistic model, (b) 4th or- 
der centered difference approximations with 12 elements on the interval [0, ir] 
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(9(7 4 ,a 2 ) holistic model is more accurate for < a < 70 but the improve- 
ment is not as pronounced as it is on the coarser grid of 8 elements. We 
suggest that this is because the major benefit to using the holistic models 
comes from application on coarser grids where the subgrid scale modelling is 
more significant. 



4.3 Comparison to Galerkin approximations 

Here we investigate the traditional Galerkin and non-linear Galerkin approx- 
imations ^3] for the Kuramoto-Sivashinsky equation (0) with the periodic 
and odd conditions (}30|) . We find the holistic models compare well with the 
Galerkin methods. While the Galerkin methods are of superior accuracy for 
solving the Kuramoto-Sivashinsky system with periodic boundary con- 
ditions, because of their global nature they lack the flexibility of the local 
nature of the holistic models. Although not explored here, this local nature of 
the holistic models empowers its use with physical boundary conditions [22 
other than periodic. 

Galerkin methods seek solutions in the form which is dominantly the 
superposition of m periodic, global modes: 

u(x, t) = b k (t) sin(kx) . (31) 

k=l 

The m-mode traditional Galerkin approximation is 

^ « (-Ak 4 + ak 2 ) b k - affi , 1 < k < m , (32) 

where 

ffl{bi, ■ ■ ■ , bm) = 2 £>i [bk +j + sign(A; - j)b lk ^] . (33) 
i=i 



The m-mode first iterate nonlinear Galerkin approximation is 
based upon the adiabatic approximation (|33j) for higher wavenumber modes 
k = m + 1 : 2m, namely 

^ w (-Ak 4 + ak 2 ) b k - aPl m (h, . . . , b m , <f> m+1 , . . . , 2m ) , (34) 
for 1 < k < m , where 

^• = -A^ m (^,---,^,0,...,0) , (35) 
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Figure 15: Bifurcation diagrams for (a) 3 mode, (b) 4 mode, (c) 6 mode and 
(d) 8 mode traditional Galerkin approximations on [0,7r]. 

Figure 16: Bifurcation diagrams for (a) 3 mode, (b) 4 mode, (c) 6 mode and 
(d) 8 mode first iterate nonlinear Galerkin approximations on [0, ir]. 

for m + 1 < j < 2m and /3 2 m is given by (J33|) . 

Obtain higher order nonlinear Galerkin approximations [22] through recog- 
nising time derivatives of these and even higher wave number modes. We do 
not explore these. 

Now examine the bifurcation diagrams of the two Galerkin approxima- 
tions (|3*TH33|) for < a < 70 and compare with the bifurcation diagrams of 
the holistic models on coarse grids, presented in M4.2I Figure shows the 
Bifurcation diagrams for the 3 mode, 4 mode, 6 mode and 8 mode traditional 
Galerkin approximations on [0, 7r]. See that at least 4 modes are needed to 
qualitatively reproduce the behaviour of the stable bimodal branch. Compare 
the 0(y 5 , a 2 ) holistic model with 6 elements from Figure ITUk. to the 6 mode 
traditional Galerkin approximation and observe the O^ 5 , a 2 ) holistic model 
qualitatively models most steady state dynamics that are reproduced by the 
6 mode traditional Galerkin approximation. Neither the (9(7 5 ,a 2 ) holis- 
tic model nor the 6 mode traditional Galerkin approximation qualitatively 
reproduce the correct behaviour of the unstable quadrimodal branch. Sim- 
ilarly the (9 (7 s , a 2 ) holistic model with 8 elements from Figure |§k and the 
8 mode traditional Galerkin approximation qualitatively model most steady 
state dynamics. However, the 8 mode traditional Galerkin approximation is 
more accurate. 

Figure EI] shows the bifurcation diagrams for the 3 mode, 4 mode, 6 mode 
and 8 mode first iterate nonlinear Galerkin approximations ((HH) on [0, it]. See 
impressive accuracy for the low mode first iterate nonlinear Galerkin approx- 
imations. The 6 mode nonlinear Galerkin approximation reproduces all of 
the steady state dynamics for the range < a < 70. There is no dis- 
cernible difference between the bifurcation diagram of the 8 mode nonlinear 
Galerkin approximation and the accurate bifurcation diagram for this range 
of a. Table E] lists the values of nonlinearity parameter a at which bifurcation 
points occur for the coarse grid holistic models and the Galerkin approxima- 
tions [T3]. The low mode first iterate nonlinear Galerkin approximations are 
impressively accurate. 

This evidence suggests that the holistic models are competitive with tra- 
ditional Galerkin approximations, but that nonlinear Galerkin models are 
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Table 2: a values at which bifurcation points occur for the various coarse 
grid holistic models and low mode Galerkin approximations 
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significantly better. However, recall that the holistic models are based upon 
analysis of local dynamics and thus we expect them to be more flexibly useful 
in applications than the global methods of these Galerkin approximations. 

4.4 Coarse grids allow large time steps 

A major benefit of accurate models on coarse grids is that larger time steps 
are possible while maintaining numerical stability. ^4. 21 shows the remarkable 
accuracy of the 0{^j 5 , a 2 ) holistic model (fl"9|) on a coarse grid of 8 elements. 
Here we investigate the maximum stable time step for explicit Runge-Kutta 
time integration on various holistic models — implicit integration schemes are 
not considered. 

In particular we compare approximations of similar accuracy but differ- 
ent grid resolutions to demonstrate the superior performance of the holistic 
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Figure 17: Bifurcation diagrams of (a) O (7 s , a 2 ) holistic model, with 8 el- 
ements on [0, 7r], and (b) 2nd order centered difference approximation with 
16 grid points on [0, 7r]. Accurate bifurcation diagram is shown in grey. 

Table 3: Approximate maximum time steps for stability of 4th order Runge- 
Kutta scheme. 
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models. For example, Figure El compares the bifurcations diagrams of the 
(9(7 5 , a 2 ) holistic model with 8 elements on [0, ir] and the 2nd order centered 
difference approximation (j2~Uj) with 16 grid points on [0,7r]. The accurate bi- 
furcation diagram is shown in grey. See that the O (7 s , a 2 ) holistic model 
on the coarse grid is of similar accuracy to the 2nd order centered difference 
approximation on the more refined grid. Thus a reasonable comparison of 
computability is made using these two schemes. 

Numerical experiments used the 4th order Runge-Kutta scheme to esti- 
mate the maximum stable time step for different holistic models and centered 
difference approximations at various values of nonlinearity parameter a. Ta- 
ble El lists the approximate maximum time steps that maintain numerical 
stability along both the negative unimodal branch at a = 10 , and the neg- 
ative bimodal branch at a = 20 and a = 30 . For the 0(y 5 ,a 2 ) holistic 
model with 8 elements, the maximum time step maintaining numerical star- 
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bility is approximately 10 times larger than the corresponding time step for 
the 2nd order centered difference approximation with 16 grid points. The 
0{^j 5 , a 2 ) holistic model requires approximately 3 times the number of float- 
ing point operations per grid value at each time step compared to the 2nd or- 
der centered difference approximation. However, on a coarse grid of 16 points 
the 2nd order centered difference approximation must be applied at twice as 
many grid points. Thus the 0(7 5 , a 2 ) holistic model can be integrated an or- 
der of magnitude faster than the 2nd order centered difference approximation 
while maintaining similar accuracy. 

Note: Table H3 shows that the higher order terms in the nonlinearity a, 
generated by the holistic method, do not reduce numerical stability. Wider 
stencil holistic approximations reduce the maximum stable time step some- 
what, but so do the wider stencil conventional centered difference approxi- 
mations. Thus, bear in mind that we need to balance the accuracy gained by 
using higher order approximation in 7, that is, wider stencil approximations, 
with the reduction in numerical stability and the increase in computation 
per grid value. 

5 Holistic models are accurate for time 
dependent phenomena 

The Kuramoto-Sivashinsky equation (JTJ) has rich dynamics [TBJ EH EE GHE El 
IT2"1 12] • Having established the excellent performance of the holistic models 
in reproducing the steady states of the Kuramoto-Sivashinsky system in 
Section EJ we now investigate the holistic models performance at reproducing 
time dependent phenomena. The Kuramoto-Sivashinsky system exhibits 
complex time dependent behaviour such as limit cycles, period doubling and 
spatio-temporal chaos. This provides us with an example to explore the 
holistic approach to modelling time dependent phenomena with relatively 
coarse discretisations. 

We restrict attention to 2n periodic solutions, 

u(x, t) = u(x + 2tt, t) . (36) 

Initially we restrict further to solutions with odd symmetry, as in the pre- 
vious section, which exhibit, see Figure |HJ Hopf bifurcations to limit cycle 
solutions, and subsequent period doubling bifurcations apparently leading to 
low-dimensional chaos E3 ED] • In §5.11 we examine the dynamics of the 
holistic models on coarse grids through the eigenvalues of the models near 
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Figure 18: The four largest (least negative) eigenvalues along the stable 
bimodal branch for the (a) 0(^f 3 ,a 2 ), (b) 0(7 4 ,a 2 ), (c) 0(j 5 ,a 2 ) holistic 
models shown in green for 8 elements on [0, 7r]. The accurate eigenvalues are 
shown in blue. 

the steady states. For example, we see that the O^ 5 , a 2 ) holistic model re- 
produces much of the eigenvalue information for < a < 70 on a coarse grid 
of 8 elements. In £ 15. 21 we explore the bifurcation diagrams near the first Hopf 
bifurcation and capture the stable limit cycles and period doubling sequence. 
The holistic models more accurately model the dynamics than centered dif- 
ference approximations of equal stencil width. Subsequently we just require 
spatial periodicity whence stable travelling wave appear followed by, at higher 
values of nonlinearity parameter a, more complex spatio-temporal chaos as 
investigated by Holmes, Lumley & Berkooz ^2] an d Dankowicz et al. jH]. 
In £ 15.31 we find the holistic discretisations more accurately model the ampli- 
tude and wave speed of travelling wave solutions, and predict better space 
time plots and time averaged power spectra, than corresponding the centered 
difference approximations. 

5.1 Dynamics near the steady states are reproduced 

Consider the eigenvalues of the Kuramoto-Sivashinsky system linearised 
about the the steady states and restricted to odd symmetry. Accurate mod- 
elling of the eigenvalues near the steady states is a necessary condition for the 
accurate modelling of the dynamics. We look at two views of the eigenvalues: 
first, their value on the negative bimodal branch; and second a more quali- 
tative plot of their values on the entire bifurcation diagram for nonlinearity 
parameter < a < 70 . 

Compare eigenvalues along the bimodal branch We investigate dy- 
namics near the stable negative bimodal branch. Consider the real part of the 
four largest (least negative real part) eigenvalues for low order holistic mod- 
els and compare to explicit centered difference approximations on a coarse 
grid of 8 elements on [0, ir]. Figure fTHI shows the four largest eigenvalues 
for the C(7 3 ,a 2 ) (JTSJ), C(7 4 ,a 2 ) (HE) and C(7 5 ,a 2 ) JEJ holistic models in 
green and the accurate solution in blue. 5 Recall the C(7 3 ,a 2 ), C(7 4 ,a 2 ) 

5 As in Section^] the accurate reference for solutions is found using a 6th order centered 
difference approximation with 48 grid points on [0, n]. 
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Figure 19: The four largest eigenvalues along the stable bimodal branch for 
the (a) 2nd order, (b) 4th order, (c) 6th order centered difference approxima- 
tions shown in magenta for 8 grid points on [0, it]. The accurate eigenvalues 
are shown in blue. 

Figure 20: Bifurcation diagram of the (9(7 5 ,a 2 ) holistic model with 8 ele- 
ments and odd symmetry on [0, 7r], depicting the real parts of the 8 largest 
(least negative) eigenvalues colour coded according to the colour bar shown. 

and (9(7 5 ,a 2 ) holistic models have 5 point, 7 point and 9 point stencils, re- 
spectively. Figure ITHh. shows the four largest eigenvalues for the C(7 5 ,a 2 ) 
holistic model closely matches the accurate solution over this range of non- 
linearity parameter a. 

Similarly, Figure shows the four largest eigenvalues for the 2nd or- 
der (|2()jl . 4th order (|2*T|) and 6th order (|22jl centered difference approximations 
in magenta on the same coarse grid. The centered difference approximations 
shown here are of equal stencil width to the corresponding holistic models 
in Figure EH Figure IT9*k . shows the 2nd order centered difference barely 
approximates the behaviour of the stable bimodal branch for a < 20 . Even 
the 6th order centered difference approximation, Figure H5h . is inferior to 
the 0(y 4 ,a 2 ^ holistic model for a > 30 . This is despite the 6th order cen- 
tered difference model having a wider stencil of 9 points compared to the 
7 point stencil of the (9(7 4 ,a 2 ) holistic model. Figures ITH1 and IT9*1 show the 
low order holistic models are superior to the corresponding centered differ- 
ence approximations for reproducing the dynamics near the stable bimodal 
branch. 

Compare eigenvalues across the bifurcation diagram Here we ex- 
plore a new view of the earlier bifurcation diagrams that additionally depicts 
the real part of the 8 largest (least negative) eigenvalues by colour. Compare 
the eigenvalues of the 0(j 5 ,a 2 ) |TH|) holistic model, see Figure E3 on the 
coarse grid of 8 elements on [0, it] to accurate ones for the Kuramoto-Sivashin- 
sky system, see Figure |^ over the nonlinearity parameter < a < 70 . The 
magnitude of the real part of the eigenvalues is colour coded according to the 

Figure 21: Bifurcation diagram of the accurate Kuramoto-Sivashinsky sys- 
tem, depicting the real parts of the 8 largest (least negative) eigenvalues, 
colour coded according to the colour bar shown. 
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Figure 22: Bifurcation diagrams near the first Hopf bifurcation for 
(a) C(7 3 ,a 2 ), (b) 0{^ 4 ,a 2 ), (c) 0(^f 5 ,a 2 ) holistic models with 8 elements 
on [0, 7r] and (d) an accurate bifurcation diagram. Stable limit cycles are 
shown in light blue and unstable limit cycles are shown in orange. 

colour bar shown on the right of the bifurcation diagram; the least negative 
eigenvalues are plotted above the more negative to give a small band of colour 
for each branch of steady states at each parameter value. Similarly to the 
bifurcation diagrams shown in Section 0J the open squares denote bifurcation 
points and the black squares denote Hopf bifurcations. Figure EDI when com- 
pared to Figure 121} shows that in addition to reproducing the stability of the 
accurate Kuramoto-Sivashinsky system for < a < 70 as discussed in §4.2^ 
the 0(7 5 , a 2 ) holistic model reproduces well the eigenvalues for most of this 
range of nonlinearity parameter a. This accurate modelling of the eigenval- 
ues is evidence of accurate modelling of theKuramoto-Sivashinsky dynamics, 
at least near the steady states. 

5.2 Extend the Hopf bifurcations 

Hopf bifurcations give rise to time periodic solutions (limit cycles). We ex- 
plore the predictions of the various models to see how well they capture these 
strongly time dependent phenomena. 

Here we investigate the bifurcation diagrams obtained by extending the 
first Hopf bifurcation, at a = 30.345 , on the positive bimodal branch and the 
period doubling sequence that ensues. We compare the bifurcation diagrams 
of low order holistic models to explicit centered difference models on a coarse 
grid of 8 elements on [0, ir] to the accurate bifurcation diagram of the Kura- 
moto-Sivashinsky system. Trajectories in the period doubling sequence are 
reported and compared by MacKenzie ^H]- As before, the holistic models 
outperform the corresponding centered difference approximations. 

Investigate the first Hopf bifurcation We now investigate the holistic 
models near the first Hopf bifurcation on the positive bimodal branch, la- 
beled HB 1; with a coarse grid of 8 elements on [0,7r]. Figure l2"2l shows the 
bifurcation diagrams of the low order holistic models and the accurate bifur- 
cation diagram near the first Hopf bifurcation. The stable limit cycles (light 
blue) that continue from this bifurcation point undertake a period doubling 
sequence commencing at a point labeled PD (yellow square). The pair of 
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Figure 23: Bifurcation diagrams near the first Hopf bifurcation for (a) 4th or- 
der, (b) 6th order centered difference approximations with 8 grid points on 
[0, 7r]. Stable limit cycles are shown in light blue and unstable limit cycles 
are shown in orange. 

unstable limit cycles born at PD give rise to the period doubling sequence 
leading to chaos. 

The accurate bifurcation diagram shown in Figure l2"2"H. is produced using 
a 6th order centered difference approximation with 24 grid points on [0, 7r]. 
The accurate bifurcation diagram shown is identical to the bifurcation dia- 
gram for the same range of a produced by Jolly et al. j!4j . Figure l2"2"k. shows 
that even the lowest order O (7 s , a 2 ) (115)) holistic model reproduces the the 
first Hopf bifurcation and finds the period doubling point on this coarse grid 
of 8 elements 6 In comparison, the corresponding 2nd order centered difference 
approximation does not even have the first Hopf bifurcation, see Figure IT2"b. 

Figurel2*2"b.c. show that higher order holistic models accurately model the 
first Hopf bifurcation and the resulting stable and unstable limit cycles. The 
accuracy of the C(7 5 , a 2 ) f|T9|) holistic model for reproducing these periodic 
solutions of the Kuramoto-Sivashinsky system is remarkable on this coarse 
grid. Figure 12*31 shows the corresponding bifurcation diagrams for the 4th or- 
der and 6th order centered difference approximations with 8 grid points on 
[0, 7r]. Compare Figure E2] and Figure 12*21 to see that the 6th order centered 
difference approximation which has a nine point stencil does not perform as 
well as the (9(7 4 ,a 2 ) (fTH|) holistic model which has a 7 point stencil. Fig- 
ure E2b,c, show that higher order holistic models more accurately model the 
first Hopf bifurcation and the resulting stable and unstable limit cycles. The 
accuracy of the C(7 5 ,a 2 ) (fE?j) holistic model for reproducing these periodic 
solutions of the Kuramoto-Sivashinsky system is remarkable on this coarse 
grid. Table |U shows the parameter values a for the Hopf bifurcations, HBi 
and the initial period doubling point PD. See that both the £?(7 4 , a 2 ) and 
(9 (7 s , a 2 ) holistic models are more accurate than the 4th order and 6th order 
centered difference approximations in reproducing the first Hopf bifurcation 
and the resulting period doubling point. 



6 Figure l2*2*k displays the bifurcation diagram for 25 < a < 32 compared to 30 < a < 37 
for the other diagrams. Since the first Hopf bifurcation for the 0(j 3 ,a 2 ) holistic model 
occurs at a = 25.595 the bifurcation diagram is shifted to contain the important dynamics. 
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Table 4: Nonlinearity parameter a values for the first Hopf bifurcation 
point HBi and resulting period doubling point PD. 
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Figure 24: a = 5 : wave-like solutions at t = 0, 0.2, 0.4, 0.6, 0.8, 1 for the 
(9 (7 s , a 2 ) holistic model shown in green and the 2nd order centered difference 
approximation in magenta on a coarse grids of 8 elements on [0,27r]. The 
accurate solution is shown in blue. 

5.3 Dynamics of periodic patterns without odd 
symmetry 

Consider the Kuramoto-Sivashinsky system (0) with solutions that are spa- 
tially periodic (|36|). We remove the requirement for odd symmetry. Conse- 
quently, we now explore travelling wave like solutions at low nonlinearity a. 
Also, we investigate the spatio-temporal chaos that occurs at higher a. 

Good performance for holistic models at low a Consider the holistic 
models of the Kuramoto-Sivashinsky system (JIJ and (}3"rjj) for nonlinearity 
parameter a = 5 and a = 10 on coarse grids of 8 elements on [0, 2tv]. Fig- 
ure l2"4l shows solutions obtained from the lowest order (9 (7 s , a 2 ) (I15j) holistic 
model for a = 5 in green, the accurate solution in blue and the corresponding 
2nd order centered difference approximation (|2()|) with 8 points on [0, 2tt], in 
magenta. The solutions are shown at time slices t = 0,0.2,0.4,0.6,0.8,1, 
starting from the half-wave initial condition of u(x, 0) = | sin(a;/2)| . See the 
(9 (7 s , a 2 ) holistic model is superior to the 2nd order centered difference ap- 
proximation on this coarse grid. In particular, the amplitude of the evolving 
wave-like solution and the wave speed are more accurately reproduced by the 
(9 (7 s , a 2 ) holistic model for a = 5 . 
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Figure 25: a = 10 : wave-like solutions at t — 0, 0.2, 0.4, 0.6, 0.8, 1 for the 
C(7 3 ,a 2 ), C(7 4 ,a 2 ) and C(7 5 ,a 2 ) holistic models shown in green, light 
green and light blue respectively and the 6th order centered difference ap- 
proximation shown in red on a coarse grids of 8 elements on [0,27r]. The 
accurate solution is shown in blue. 

Similarly, Figure ESI shows the same time slices for larger a = 10 . The 
(9 (7 s , a 2 ) (|T5Jl holistic model is shown in green, the C(7 4 , a 2 ) (fTSj) model is 
shown in light green and the (9(7 5 ,a 2 ) (|19|) holistic model in light blue for 
this coarse grid of 8 elements on [0,27r]. For this a the 2nd order (|2()j) and 
4th order ()21j1 centered difference approximations do not generate a wave-like 
solution at all. However, the 6th order centered difference approximation (|2"2"jl 
does produce the travelling wave-like solution shown in red. The (D(y 3 ,a 2 ) 
holistic model (green) is the least accurate on this coarse grid but it does 
reproduce a stable solution on this coarse grid for only a 5 point stencil 
approximation. The (9(7 4 ,a 2 ) holistic model (light green) more accurately 
models the amplitude of the solution compared to the 6th order centered dif- 
ference approximation despite having a smaller stencil width. The 0(^j 5 , a 2 ) 
holistic model is the most accurate at reproducing both the amplitude and 
the wave speed of the stable wave-like solution for a = 10 on this coarse grid 
of 8 elements. 

Good performance for more complex behaviour For higher values 
of nonlinearity parameter a for which the Kuramoto-Sivashinsky system ex- 
hibits more complex behaviour, including spatio-temporal chaos, it is more 
useful to compare time averaged power spectra rather than particular trav- 
elling waves. Here we investigate the performance of the holistic models on 
coarse grids for a = 20 and 50 using the example of the 0(7 5 , a 2 ) (|19|) holis- 
tic model on relatively coarse grids, and we compare it with the 6th order 
centered difference approximation which is of equal stencil width. The focus 
here is to show the improved performance of the holistic models for ranges of 
parameter a that contain more complex time dependent behaviour. We also 
expect a corresponding improvement for the other holistic models but this is 
not investigated here. Further, we also compare the O^ 5 , a 2 ) holistic model 
on coarse grids to the 2nd order centered difference approximations of similar 
accuracy. We find the 0(y 5 ,a 2 ) holistic model, but with approximately 1/3 
of the grid points, has comparable accuracy to 2nd order centered difference 
appr oximat ions . 

Figure EH1 shows space time plots of (a) the 0(7 5 ,a 2 ) holistic model 
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Figure 26: a = 20 : space time plots for (a) the (9(7 5 ,a 2 ) holistic model 
with 12 elements on [0, 27r], (b) 6th order centered difference approximation 
with 12 grid points on [0, 2ti} and (c) the accurate solution 

Figure 27: a = 20 : time averaged power spectra for the 0(7 5 ,a 2 ) holistic 
model with 12 elements on [0, 2tt] shown in light blue, and the 6th order 
centered difference approximation in red for (a) 12 grid points on [0, 2-k] and 
(b) 16 grid points on [0, 2ir}. The accurate power spectrum is shown in blue. 

with 12 elements on [0, 2tt], (b) the 6th order centered difference approxima- 
tion with 12 grid points on [0, 2n} and (c) the accurate solution. 7 See the 
0(y 5 , a 2 ) holistic model reproduces much of the complex structure of the ac- 
curate solution for a = 20 with 12 elements. Figure l2*^b. shows the 6th order 
centered difference approximation incorrectly finds a periodic solution after 
approximately t = 0.2 . 

Since the Kuramoto-Sivashinsky system at nonlinearity parameter a = 
20 exhibits more complex time dependent behaviour than simple limit cy- 
cles, we compare time averaged power spectra, denoted here by S(k) for 
wavenumber k. Figure l2*Tk . shows a log-log plot of the time average power 
spectra of the (9 (7 s , a 2 ) holistic model in light blue and the 6th order cen- 
tered difference approximation on a coarse grid of 12 elements on [0, 2ii\ in 
red. The accurate power spectrum is shown in blue. For this coarse grid of 
only 12 elements only 5 wavenumbers are relevant, as displayed. See that 
the (9(7 5 ,a 2 ) holistic model is superior to the 6th order centered difference 
approximation on this coarse grid of 12 elements. Figure l2*7b. compares the 
time average power spectrum of the 0(j 5 ,a 2 ) holistic model in light blue 
with 12 elements and the 6th order centered difference approximation with 
16 grid points. The (9(7 5 ,« 2 ) holistic model achieves similar accuracy on a 
coarser grid. 

7 The accurate solutions plotted in this section are computed using a 6th order centered 
difference approximation and 256 grid points on the interval [0, 2tt]. This is sufficient grid 
resolution to capture the important dynamics of the Kuramoto-Sivashinsky system for 
the values of a investigated here. 

Figure 28: a = 20 : time averaged power spectra for for the O^ 5 , a 2 ) holistic 
model with 12 elements on [0, 2n] shown in light blue, and the 2nd order 
centered difference approximation in magenta for (a) 24 grid points on [0, 2tt] 
and (b) 36 grid points on [0, 27r]. The accurate spectrum is shown in blue. 
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Figure 29: a = 50 : space time plots for (a) the 0(y 5 ,a 2 ^ holistic model 
with 24 elements on [0, 27r], (b) 6th order centered difference approximation 
with 24 grid points on [0, 2tt] and (c) the accurate solution. 

Figure 30: a = 50 : time averaged power spectra for the 0[j 5 ,a 2 ) holistic 
model with 24 elements on [0, 2tt] shown in light blue, and the 6th order 
centered difference approximation in red for (a) 24 grid points on [0, 2-k] and 
(b) 32 grid points on [0, 2%]. The accurate spectrum is shown in blue. 

The power spectra of the 0(7 5 ,a 2 ) holistic model on a coarse grid of 
12 elements and the 2nd order centered difference approximation on the more 
refined grids of 24 and 36 points are shown in Figures l2*Tk.b respectively. See 
a refined grid of 36 points is needed achieve similar accuracy to the 0(y 5 , a 2 ) 
holistic model on a coarse grid of 12 elements on [0, 2ir]. That is, through 
its subgrid scale modeling, the holistic model achieves similar accuracy with 
one-third the dimensionality. 

For nonlinearity parameter a = 50 the Kuramoto-Sivashinsky system 
exhibits even more complex behaviour, see the space time plots in Figure l29l 
The 0(^y 5 ,a 2 ) holistic model more accurately reproduces the Kuramoto- 
Sivashinsky system than the 6th order centered difference approximation on 
this coarse grid of 24 elements. On this coarse grid the 6th order centered 
difference approximation shown in Figure l2*9"b . exhibits a periodic solution 
after time t ^ 0.1 which does not match the irregular behaviour seen in the 
accurate solution and the C(7 5 ,a 2 ) holistic model. 

We again examine time averaged power spectra to further investigate the 
performance of the 0{^j 5 , a 2 ) holistic model at this relatively large parameter 
value of a = 50. Figure EH compares the time averaged power spectrum of the 
0(y 5 , a 2 ) holistic model in light blue on a coarse grid of 24 elements on [0, 2ix] 
to the 6th order centered difference approximation in red, for (a) 24 grid 
points and (b) 32 grid points on [0, 2tt]. The 6th order centered difference 
approximation with 32 grid points has similar accuracy to the C(7 5 ,a 2 ) 
holistic model on a coarse grid of just 24 elements for a = 50 . 

This investigation of the O (7 s , a 2 ) holistic model on coarse grids for 
a = 20 and 50 shows it reproduces similar accuracy to the 2nd order centered 
difference approximation on a coarse grid of approximately 1/3 the resolution, 
and similar accuracy to the 6th order centered difference approximation on 
grids of approximately 3/4 the resolution. MacKenzie ^H] reports that even 
at a — 200 the holistic model qualitatively well captures the dynamics of the 
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Kuramoto-Sivashinsky pde. This increased accuracy on coarse grids allows 
larger time steps for explicit time integration schemes, as discussed in £ 14.41 

6 Conclusion 

Holistic discretisation [23] is straightforwardly extended to fourth order dissi- 
pative PDEs through the example of the Kuramoto-Sivashinsky equation ^H] • 
We divide the domain into elements by introducing artificial internal bound- 
ary conditions (§2J) which isolate the elements when 7 = but when 7 = 1 
they fully couple the elements to recover the Kuramoto-Sivashinsky equa- 
tion. Then centre manifold theory supports the discretisation, see £ 12.21 The 
holistic models listed in §Hlhave a dual justification ( §3.3|) : not only are they 
supported by centre manifold theory for finite element size h, the ibcs are 
specially crafted |2E] so the models are also consistent with the Kuramoto- 
Sivashinsky equation as the grid spacing h — > . 

No formal error bounds currently exist for the holistic method; the diffi- 
culty is that the models are based at 7 = but are evaluated at finite 7 = 1. 
Instead we present a detailed numerical investigation of the holistic models 
of the steady states (Section HJ) and time dependent solutions (Section 0) of 
the Kuramoto-Sivashinsky on coarse grids. 

We compared, in H4A\ the accuracy of different approximations in pre- 
dicting steady states on different grid resolutions. The holistic O (7 s , a 2 ) 
approximation on a grid of 8 elements has similar accuracy to a 2nd order 
centered difference approximation on a grid of 16 points. Consequently the 
holistic model allows a maximum time step which is an order of magnitude 
longer than that of the explicit centered difference approximation of similar 
accuracy, while maintaining numerical stability. The accuracy of the holistic 
approximations to the Kuramoto-Sivashinsky equation on coarse grids and 
subsequent improved performance justifies further application of the holistic 
method and future investigation of the approach. 

The holistic models on coarse grids also modelled well time dependent 
phenomena of the Kuramoto-Sivashinsky system. In particular, in £ 15.11 we 
saw the holistic models more accurately model the eigenvalues near the steady 
states of the first form of the Kuramoto-Sivashinsky system compared to ex- 
plicit centered difference approximations of equal stencil widths. The coarse 
grid holistic models also more accurately model the first Hopf bifurcation 
and the resulting period doubling sequence, see £ 15.21 Further, in comparison 
with explicit centered difference models, in §5. 31 we saw good performance for 
higher values of the nonlinearity parameter a and more accurate predictions 
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of time averaged power spectra: the 0(^j 5 , a 2 ) holistic model achieves similar 
accuracy to the 2nd order and 6th order centered difference approximations 
on approximately 1/3 and 3/4 of the grid resolutions respectively. 

This good performance of the holistic models for accurately reproducing 
both the steady states and the time dependent phenomena of the Kuramoto- 
Sivashinsky system is good evidence that the holistic approach is a powerful 
method for discretising dissipative pdes on coarse grids. 

A Computer algebra derives the 
discretisation 

1 write "The holsitic discretisation of a supplied PDE. "$ 

2 write "based on holistic PDE version beta.l, 20 Dec 2001"$ 
3 

4 stenwidth:=9; 

5 epsilo:=2; 
6 

7 % improve printing 

8 linelength 72$ 

9 on div; off allfac; on revpri; 
10 factor gamma, h; 

11 

12 °/ make function of xi=(x-x_j)/h 

13 depend xi,x; 

14 let df (xi,x)=>l/h; 
15 

16 % get parameters of the PDE 

17 write "The parameters of the PDE follow"$ 

18 operator uu; 

19 depend uu,x; 

20 dissipate : =df (uu,x, 4) ; 

21 discof:=-4; 

22 °/ set neglected orders of errors 

23 if stenwidth=3 then let gamma"2=>0 

24 else if stenwidth=5 then let gamma~3=>0 

25 else if stenwidth=7 then let gamma~4=>0 

26 else if stenwidth=9 then let gamma"5=>0 

27 else if stenwidth=ll then let gamma"6=>0 

28 else if stenwidth=13 then let gamma~7=>0 

29 else if stenwidth=15 then let gamma~8=>0 

30 else if stenwidth=17 then let gamma"9=>0 
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31 else if stenwidth=27 then let gamma" 15=>0 

32 else let gamma=>0$ 

33 if epsilo=decreasing then epsilon:=gamma 

34 else if epsilo=2 then let a~2=>0 

35 else if epsilo=3 then let a~3=>0 

36 else if epsilo=4 then let a"4=>0 

37 else if epsilo=5 then let a~5=>0 

38 else let a=>0$ 
39 

40 %decreasing:=l$ 

41 operator u; 
42 

43 % solvability condition 

44 operator solg; linear solg; 

45 let { solg(xi-~p,xi)=>(l+(-l)-p)/(p+2)/(p+l) 

46 , solg(xi,xi)=>0, solg(l ,xi)=>l }; 
47 

48 % define solving operator depending upon the dissipation 

49 operator solv; linear solv; 

50 if sub (uu=x~9, dissipate) =72*x~7 then begin 

51 write disorder: =2; 

52 % solves v"=RHS s.t. v(0)=0 and v(+l)=v(-l) 

53 let { solv(xi~~p,xi) => 

54 ( xi-(p+2)-(l-(-l)-p)*xi/2 )/(p+l)/(p+2) 

55 , solv(xi,xi) => (xi"3-xi)/6 

56 , solv(l.xi) => (xi~2)/2 >; 

57 end else 

58 if sub(uu=x~9,dissipate)=3024*x~5 then begin 

59 write disorder: =4; 

60 % solves v""=RHS s.t. v(0)=v(+l)=v(-l)=0 and v(+2)=v(-2) 

61 let { solv(xi~~p,xi) => 

62 ( xi-(p+4)-(l+(-l)~p)/2*xi-2 

63 -(l-(-l)-p)/6*((2-(p+3)-l)*xi-3+(4-2-(p+3))*xi) 

64 )/(p+l)/(p+2)/(p+3)/(p+4) 

65 , solv(xi.xi) => (xi"5- (15*xi"3-12*xi) /3) /120 

66 , solv(l.xi) => (xi"4-xi"2)/24 }; 

67 end else 

68 if sub(uu=x~9,dissipate)=60480*x~3 then begin 

69 write disorder: =6; 

70 % solves v"""=RHS s.t. v(0)=v(+l)=v(-l)=v(+2)=v(-2)=0 and v(+3)=v(- 

71 let { solv(xi"~p,xi) => 

72 ( xi"(p+6)+(l+(-l)"p)/6*((2-(p+4)-4)*xi"2+(l-2"(p+4))*xi-4) 

73 +(!-(-!) "p) /240* (4* (-45+9*2" (p+6) -3" (p+6) ) *xi 
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74 +5* (13-8*2" (p+6) +3" (p+6) ) *xi~3 

75 +(-5+4*2~(p+6)-3~(p+6))*xi"5) 

76 ) / (p+1) / (p+2) / (p+3) / (p+4) / (p+5) / (p+6) 

77 , solv(xi.xi) => (xi~7-36*xi+49*xi~3-14*xi~5) /5040 

78 , solv(l.xi) => (xi~6+4*xi~2-5*xi~4)/720 >; 



79 end; 

80 if disorder>stenwidth then 

81 write "*** Warning: the stencil width is too small ***"; 
82 

83 

84 % parametrise with evolving uu(j) 

85 depend u,t; 

86 let df (u(~k) ,t)=>sub(j=k,gj) ; 
87 

88 % linear solution in jth element 

89 write "Start with the linear approximation" $ 

90 uu:=u(j)+udash; 

91 udash:=0; 

92 gj:=0; 

93 % iterative refinement to specified error 

94 write "Iterate to make residuals negligible "$ 

95 iteration :=0$ 

96 let {u(j+4)=>-u(j+3) ,u(j+5)=>-u(j+2) ,u(j+6)=>-u(j+l) 

97 ,u(j+7)=>-u(j) ,u(j-l)=>-u(j) ,u(j-2)=>-u(j+l) 

98 ,u(j-3)=>-u(j+2) ,u(j-4)=>-u(j+3)>; 
99 

100 repeat begin 



101 write iteration :=iteration+l; 

102 deq: =-df (uu, t)+discof *dissipate ; °/ -a* (df (uu,x,2) ) ; 

103 °/,uu*df (uu,x) + 

104 rbc :=-(sub(xi=+l ,uu)-sub(xi=0,uu) )+gamma*(u(j+l)-u(j) ) ; 

105 lbc : =-(sub(xi=0 ,uu)-sub(xi=-l ,uu) )+gamma*(u( j )-u(j-l)) ; 

106 ok:= if (deq=0)and(rbc=0)and(lbc=0) then 1 else 0; 

107 if disorder>3 then begin 

108 rrbc : =-(sub(xi=+2 ,uu)-3*sub(xi=+l ,uu)+3*sub(xi=0,uu)-sub(xi=-l ,uu) ) 

109 +gamma~2*(u(j+2)-3*u(j+l)+3*u(j)-u(j-l)) ; 

110 llbc : =- ( sub (xi=+l , uu) -3*sub (xi=0 , uu) +3*sub (xi=- 1 , uu) -sub (xi=-2 , uu) ) 

111 +gamma~2*(u(j+l)-3*u(j)+3*u(j-l)-u(j-2)) ; 

112 ok:=if ok and(rrbc=0) and(llbc=0) then 1 else 0; 

113 if disorder>5 then begin 

1 14 rrrbc : =- ( sub (xi=+3 , uu) -5*sub (xi=+2 , uu) +10*sub (xi=+l , uu) 

115 -10*sub (xi=0 , uu) +5*sub (xi=-l , uu) -sub (xi=-2 , uu) ) 

116 +gamma~3*(u(j+3)-5*u(j+2)+10*u(j+l)-10*u(j)+5*u(j-l)-u(j-2)) ; 
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117 lllbc:=-(sub(xi=+2,uu)-5*sub(xi=+l,uu)+10*sub(xi=+0,uu) 

118 -10*sub (xi=-l , uu) +5*sub (xi=-2 , uu) -sub (xi=-3 , uu) ) 

119 +gamma~3*(u(j+2)-5*u(j+l)+10*u(j)-10*u(j-l)+5*u(j-2)-u(j-3)) ; 

120 ok:=if ok and(rrrbc=0)and(lllbc=0) then 1 else 0; 

121 end; 

122 end; 

123 gd:=solg(deq,xi) +discof /h~disorder* 

124 (if disorder=2 then (rbc-lbc) else 

125 if disorder=4 then (rrbc-llbc) else 

126 if disorder=6 then (rrrbc-lllbc) ) ; 

127 gj:=gj+gd; 

128 udash: =udash+h~disorder*solv(-deq+gd,xi) /discof 

129 +(if disorder=2 then xi/2* (rbc+lbc) else 

130 if disorder=4 then xi/2* (rbc+lbc)+xi~2/2* (rbc-lbc) 

131 -(xi-xi~3)/12*(rrbc+llbc) else 

132 if disorder=6 then xi/2* (rbc+lbc) +xi~2/2* (rbc-lbc) 

133 -(xi-xi"3)/12*(rrbc+llbc) -(xi~2-xi~4)/24*(rrbc-llbc) 

134 +(4*xi-5*xi~3+xi~5)/240*(rrrbc+lllbc) ); 

135 showtime; 



136 end until ok or (iteration>25) ; 
137 

138 write deq:=deq; 

139 write rbc:=rbc; 

140 write lbc:=lbc; 

141 write rrbc:=rbc; 

142 write llbc:=lbc; 

143 write rrrbc:=rbc; 

144 write lllbc:=lbc; 
145 

146 end; 
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